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The accurate and efficient description of strongly correlated systems remains an important chal¬ 
lenge for computational methods. Doubly occupied configuration interaction (DOCI), in which all 
electrons are paired and no correlations which break these pairs are permitted, can in many cases 
provide an accurate account of strong correlations, albeit at combinatorial computational cost. Re¬ 
cently, there has been significant interest in a method we refer to as pair coupled cluster doubles 
(pCCD), a variant of coupled cluster doubles in which the electrons are paired. This is simply 
because pCCD provides energies nearly identical to those of DOCI, but at mean-field computa¬ 
tional cost (disregarding the cost of the two-electron integral transformation). Here, we introduce 
the more complete pair extended coupled cluster doubles (pECCD) approach which, like pCCD, 
has mean-field cost and reproduces DOCI energetically. We show that unlike pCCD, pECCD also 
reproduces the DOCI wave function with high accuracy. Moreoever, pECCD yields sensible albeit 
inexact results even for attractive interactions where pCCD breaks down. 
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I. INTRODUCTION 


While the desciption of the ground state of weakly correlated systems is by now fairly routine, the same cannot be 
said for strongly correlated problems. Because coupled cluster theorjC^^^ offers exceptional accuracy for the description 
of weak dynamic correlations, we would like to use some variant of coupled cluster theory as well for the strong static 
correlations; in this way one could seamlessly merge the two ideas to provide a more powerful technique which should 
be accurate for both regimes. Unfortunately, the construction of coupled cluster techniques for strongly correlated 
systems is a work in progress. 

A particularly interesting recent development is the notion of pair coupled cluster doubles (pCCD),®^ which takes 
the simple coupled cluster wave function and makes the dramatic simplification that the only allowed excitations are of 
a paired form in which two electrons are removed from the same spatial orbital and placed in some other spatial orbital. 
In other words, pCCD is coupled cluster doubles restricted to include only seniority zero determinants, where the 
seniority of a determinant is the number of singly occupied spatial orbitals. This restriction greatly decreases the cost 
of the coupled cluster calculation to mean field or 0{M^) if one ignores the two-electron integral transformation, but 
rather paradoxically, despite this simplification the pCCD wave function yields results very close to the doubly occupied 
configuration interaction (DOCI),Q2Hini-vvhich includes such pair excitations to all excitation levels. This DOCI method 
is not new, and includes many powerful geminal wave functions, including the antisymmetrized geminal power (APG), 
the antisymmetric product of strongly orthogonal geminals (APSG), and many others. For many problems of interest, 
DOCI is able to describe the basics of the strong correlations, and to the extent that pCCD reproduces DOCI, so too 
does pCCD. 

We note, however, that the coincidence between pCCD and DOCI is not entirely universal. For the attractive 
pairing Hamiltonian 


iJ — ^ Cp 


(4, + 4, Op^) 


aq, 


( 1 ) 


we observe that as the pairing strength G becomes larger and the mean-field solution develops an instability toward 
a number symmetry broken Hartree-Fock-Bogoliubov state, pCCD breaks down dramatically, overcorrelating wildly 
before eventually returning complex energies.l^ This suggests that pCCD may not be able to describe the kinds of 
strong correlations needed to model superconductivity, for example. 

Moreover, even when the pCCD energy is accurate, the wave functions may be less so. In coupled cluster theory, 
the Hamiltonian is similarity transformed, yielding a non-Hermitian effective Hamiltonian H = exp(—T) H exp(T). 
Because H is non-Hermitian, it has different left- and right-hand eigenvectors: 

H\0) = UjO), (2a) 

{Q\{1 + Z)H = E{0\{1 +Z) (2b) 


where Z creates excitations to the left and |0) is the mean-field reference. This structure in turn translates to different 
left- and right-hand wave functions for the original Hamiltonian, which in pCCD are taken to be 


(-^pCCdI = (0|(I-I-2') e 
li^pcco) = e^|0). 


(3a) 

(3b) 


where T and Z respectively create pair excitations and pair de-excitations when acting to the right. In Ref. [5] we and 
others showed that the overlap 


^ = (/:pcCD|DOCI)(DOCI|7^pccD) 


( 4 ) 


is close to unity, but this does not separately test the pCCD left-hand and right-hand wave functions and, as we will 
see below, while the pCCD right-hand wave function is close to the DOCI state, the same is not necessarily true of 
the pCCD left-hand wave function. This can cause the pCCD density matrices to be less accurate approximations of 
the actual DOCI density matrices than we would like. 

This manuscript introduces the pair extended coupled cluster method (pECCD) which seeks to remedy these 
deficiencies. Where pCCD is the seniority zero version of coupled cluster doubles, pECCD is the seniority zero version 
of the extended coupled cluster doubles method of Arponen and BishopWhile extended coupled cluster has not 
seen a great deal of use due to its large computational cost (C>(M^°) for extended coupled cluster double^^, pECCD 
has the same mild scaling with system size displayed by pCCD, though with a rather larger prefactor. We will sketch 
the pCCD method in Sec. |ll]and provide a few results in Sec. |HI| to show the basic performance of the method before 
providing our conclusions and prospects for further development in Sec. IV We should note that the equations needed 



3 


for efficient computation of the pECCD energy and wave function amplitudes are exceedingly lengthy despite their 
low computational cost; here, we provide a simpler but computationally more demanding expression for the energy 
(from which the amplitude equations follow) and have not included the amplitude equations at all. Accordingly, the 
code needed to solve pECCD more efficiently will be made available upon request. 


II. THEORY 

The basics of the pair extended coupled cluster method are simple. We define a pair excitation operator 

T = Y,UaPlP. (5) 

ai 

and a pair de-excitation operator 

z = Y.^-^p^p- (6) 

ai 

where the pair creation operator is given by 

Pl = « ( 7 ) 

and the pair annihilation operator Pp is its adjoint. Note that these operators are nilpotent (Pp = 0) and that more 
general pairing schemes than this simple singlet pairing within the same spatial orbital are possible. 

Regardless, having defined an excitation operator T and a de-excitation operator Z, the pECCD energy is given by 


ApEccD = (0|e^e H e ^|0) 

= {{)\{l + zp\z^ + l Z^) He^\G) 

1 D 


(8a) 

(8b) 


where in the second line we have used the facts that the de-excitation operator Z annihilates the vacuum to the right 
and that the similarity transformed Hamiltonian H = e“^ i? e^ is a six-body operator that creates up to hextuple 
excitations to the right, so that no more than hextuple de-excitations (created by Z^) are needed for the case under 
consideration. The amplitudes Ua and Zai are obtained by solving 


0 = 


i^-^pECCD 

dtia 


^ApECCD 

dZai 


(9) 


We will take the physical Hamiltonian to be 


H = “L + \ 

pq a 


^{P<l\v\rs) al^ a\^ Qr, 

pqrs (77] 


( 10 ) 


where {p\h\q) and {pq\v\rs) are one- and two-electron integrals; note that the two-electron integrals here are not 
antisymmetrized. Greek letters index spin while Latin letters index spatial orbitals, with denoting occupied 

orbitals, a,b,c,... denoting virtual orbitals, and p,q,r,... denoting arbitrary orbitals. 

We can greatly simplify the derivation of the pECCD energy (and therefore amplitude equations) by replacing the 
physical Hamiltonian H with the portion which preserves seniority. For a two-body Hamiltonian, we would 

have 


^ ^ jjsn=o ^ ^sn=2 ^ ^SQ=4 

where iH5n=2 ^<50=4 determinants whose seniorities differ by two and by four, respectively. Determinants 

differing by an odd seniority have different electron numbers, while determinants differing by an even seniority greater 
than four differ by a triple excitation or higher (and thus cannot be coupled by a two-body operator). Because every 
spatial orbital in pCCD or pECCD has seniority zero {i.e. all orbitals are either doubly occupied or empty), it will 
suffice for our purposes to determine what we will call which is the part of that preserves the seniority 

of each orbital. We caution, however, that omits terms which change the seniority of individual orbitals while 
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preserving the total seniority; these terms do not contribute to closed-shell pCCD or pECCD but could contribute to 
a kind of ROHF-based generalization of pCCD or pECCD. 

We derive our expression for in the appendix and merely quote it here: 

= + ^p^ Pp + E ^p^ ^p ■ (12) 

P p^q pq p^q 

where the necessary integrals are 


hp = {p\h\p), (13a) 

Vpq = {pp\v\qq), (13b) 

Wpg = 2 {pq\v\pq) - {pq\v\qp), (13c) 

Kpq = -{pq\v\qp) (13d) 

and where the number and spin operators are given by 

Np = 'y ( ap^, (14a) 

cr 

^ E “p. (14b) 


where a is the vector of Pauli matrices. The Heisenberg-like term J2p^q ^pq ^p ' ^q '"^bl not contribute in our closed- 
shell case because closed-shell orbitals have spin zero. The number operator and pair creation and annihilation 
operators satisfy SU(2) commutation relationships 


{Pp,Pl]=Spq i^-Np), 

[Np,Pg] = -2SpgPg. 


One can verify that 


(O|^5n=0|o) = (0|Ef|0) 


(15a) 

(15b) 


(16) 


where |0) is an RHF or ROHE determinant (or, more precisely, any determinant which uses the same spatial orbitals 
for different spins). 

Given the seniority-preserving Hamiltonian, the pECCD energy can be obtained simply by computing the pECCD 
density matrices. As discussed in Ref. [HI the density matrices of seniority zero methods are sparse, and in fact the 
only non-zero elements are precisely those we need to evaluate the expectation value of We may of course 

also use density matrices to take other expectation values; in particular, they can be used in orbital optimization.^^ 
In this work, we will use orbitals optimized for pCCD, because in our experience the pCCD- and pECCD-optimized 
orbitals are virtually identical and the pCCD orbital optimization is somewhat less expensive because the pCCD 
density matrices are simpler to compute. 

To compute density matrix elements, we define a pECCD expectation value via 


(O) = (0|e^e-^C>e^|0). 


(17) 
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Given this expectation value, we find that the density matrices are given by 


{N,) = 2 

{Na)=2J2t^a 


ta ^a% I 5 


^ I ^ ^ ^ {^ia ^ai ^ja ^aj) ^ij ^ ^ ^ia ^ai H" ^ij ^ ^ ^ia ^jb i^ai ^bj ^bi ^aj) | ; 


a^b 


i^a ^b) — 4 I Sab ^ ^ ^ai ^ia ^ab ^ ^ ^ia ^jb i^ai ^bj H“ ^bi ^aj) | : 




-^a) — 4 I ^ ^ t^a ^ai ^ ^ ^ ^ ^ia ^jb {_^ai ^bj H“ ^bi ^aj) | ; 
\i^j i^j b^a 

{Pa Pi) = Zaii 

{Pa Pb) — ^ ^ ^ai ^ib 2 Sab ^ ^ ^ai ^bj ^ib ^jb^ 




{Pi Pj) — ^ ^ ^ia ^aj H“ Sij 11 2 ^ ^ tia Zai I 2 Sij ^ ^ t^a Zai Z})j , 

a \ a / a^h 

{Pi Pa) — ^ia 2 ^ ^ tia tja Z-aj 2 ^ ^ til) tia Zbi ^ ^ tib tja Zbj 2 ^ ^ ^ ^ tib tja tf^a Zaj Zbk 

j^i b bj b^a j^k 

H” 4 ^ ^ ^ ^ tia tib tja {Zai Zbj “h Zbi Zaj ) 2 ^ ^ ^ ^ tib tic tja Zbi Zcj 

ji^i b^c j^i 

+ 2 E E ^ib tic tja tka ^bj ^ck “t 2 Zaj Zl)k Zci) 

b-^c^a j^k^i 


(18a) 

(18b) 

(18c) 

(18d) 

(18e) 

(18f) 

(18g) 

(18h) 

(18i) 


where 5pq = 1 — 5pq. Then the energy is just given by 

^^pECCD = E + E ^ ^ ^ E + b ^ 

a¥^b 

+ E {pI Pj) + E P-) + + E (Pa Pb) 


ja 


where we have disregarded the Heisenberg Hamiltonian term as its expectation value vanishes. One can verify that 
the pCCD energy is properly reproduced by taking the pECCD energy and omitting terms of 0{Z^) or 0{Z^). 

While constructing the density matrices given in Eqn. 18 would appear to require 0{N^) computational time, 
the summation restrictions can be lifted by adding and removing terms which, with sufficient exertion, allows one to 
evaluate the density matrices and thus the energy in 0{N^) time after introducing intermediates. For example, we 
may write 


in terms of intermediates 


tia tib ^ai Zbj — E tia tib Zai Zbj -E tia tia Zai Za 


a^b 


— ^ii ^ij ^ij 


^ij — ^ ^ tia Zc 
a 

S^ij — ^ ^ tia ti, 


Zai Zai ■ 


(20a) 

(20b) 


(21a) 

(21b) 
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Similarly, the amplitude equations can be solved in 0{N^) operations with the appropriate definition of intermediates. 
We have checked the correctness of our 0{N^) implementation by comparison to the explicit 0{N^) result for random 
input T and Z amplitudes and integrals h, v, and w, and have verified our 0{N^) amplitude equations by comparing 
analytic and numerical derivatives of ifpECCD for random inputs. 

We should emphasize that while the pCCD and pECCD density matrices both adopt a quasi-diagonal form, the 
pECCD two-particle density matrices are much more complicated. This is simply because the pECCD left-hand wave 
function (see below) is more sophisticated; it contains excitations to all even orders and thereby has more flexibility 
in fitting DOCI than does pCCD, even though pCCD and pECCD have the same number of parameters to optimize. 


III. RESULTS 

Following Ref. El we will compare the pCCD, pECCD, and DOCI energies for a variety of systems, defining, for 
example. 


^^'pECCD = EpECCD — EdOCI) (22) 

and will also assess the quality of the pCCD and pECCD wave functions by evaluating 

S = (riDOCI) (DOCI|7^), (23) 

where for both pCCD and pECCD the right-hand wave function is 

1^) = e^|0) (24) 

while the left-hand wave functions are 

(>Cpeccd| = (0|e^e“^ 

(-CpcCDl = (0|(1 -I- Z) e~'^ 

for pECCD and pCCD, respectively. Recall that 

1 = {C\TZ) = (riDOCI) (DOCI|7^) + ^(£|DOCIfc) (DOCIfc|7^) (26) 

k 

where |DOCIfe) is the excited DOCI state; thus, we will have 5” « 1 provided that the DOCI excited states have 
minimal overlap with either the left-hand or right-hand state of pECCD or pCCD. It will also prove fruitful to look 
in more detail, however, at the individual left- and right-hand overlaps 

=AfL(£|DOCI), (27a) 

5i^=Affl(DOCI|7^), (27b) 

where the normalization constants Ml and Mu are such that the left- and right-hand states are individually normalized 
to unity so that, for example. 


(25a) 

(25b) 


Ml{n\n)=Ml{Q\e^'e^\Q) = l. (28) 

These allow us to assess separately the quality of the left- and right-hand wave functions. 

As Ref. Elniakes clear, with orbital optimization, pCCD is exact for two-electron singlets. The same is of course 
true for pECCD. Thus, both are exact for H 2 and nearly exact for LiH, with little to distinguish the two approaches 
in the latter case. It will therefore be more fruitful to focus on systems with more strongly correlated electrons. 

We begin, then, with chains of equally spaced hydrogen atoms. The strong correlations in these systems seem to 
be described reasonably well by DOCI. Figure shows that while pCCD reproduces DOCI quite well energetically, 
pECCD does so even better. More relevantly, while both pCCD and pECCD accurately describe the right-hand wave 
function (ie. we see that |DOCI) « exp(T)|0)), the left-hand wave function of pCCD is a fairly poor approximation 
to the DOCI state, while the left-hand wave function of pECCD is again very accurate. This should not be too 
surprising, as the pCCD left-hand wave function (0|(1 -I- Z) exp(—T) consists only of the reference and doubly excited 
determinants, while the pECCD left-hand wave function includes all the same higher excitations that DOCI adds. 
Indeed, that pECCD accurately reproduces the left-hand wave function of DOCI while pCCD sometimes does not 
seems to be a fairly general feature. Note that while an accurate left-hand wave function is not needed for accurate 
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FIG. 1. Differences between pCCD, pECCD, and DOCI in equally spaced hydrogen chains. Top left: Base 10 logarithm of the 
absolute value of AS, measured in Hartrees and defined in Eqn. |22| Top right: Base 10 logarithm of the absolute value of 
1 — S' as defined in Eqn. |23| Bottom left: Base 10 logarithm of the absolute value of 1 — Sn as defined in Eqn. |27a[ Bottom 
right: Base 10 logarithm of the absolute value of 1 — Sr as defined in Eqn. |27b| Kinks in the pCCD results for AE and 1 — S 
are due to changes in sign. 


energies, errors in the left-hand wave function may translate into errors for properties other than the energy. In other 
words, while both pCCD and pECCD match the DOCI energy, we would expect only pECCD to match DOCI for 
arbitrary properties. We should also perhaps emphasize the smallness of the various discrepancies; the pECCD energy 
differs from DOCI by less than 0.001 kcal/mol per electron, and the left- and right-hand states have overlaps with 
DOCI differing from I by about I0“®, implying that the coefficients of DOCI excited states in the pECCD ground 
state are less than 10“"^. For all practical purposes, pECCD reproduces DOCI exactly. 

Similar results are seen in the symmetric double dissociation of H 2 O (Fig. and in the dissociation of N 2 (Fig. |^. 
While the pCCD energy and right-hand wave function are close to those of DOCI, pECCD is closer yet; meanwhile, 
the pECCD left-hand wave function may be a much better approximation to DOCI than is the pCCD left-hand 
wave function. This is true not just in the basis of energetically optimized orbitals, but is also true with canonical 
Hartree-Fock orbitals (see, e.g., Fig. |^. 

We have mentioned that the close coincidence between pCCD and DOCI breaks down for the attractive pairing 
Hamiltonian of Eqn. [^which in the language of the SU(2) generators discussed earlier is just 

H = (29) 

pq 

Because this Hamiltonian contains only the SU(2) pairing or pseudospin genera tors, it is solved exactly by DOCI. 
Conveniently, however, a more compact exact solution was found by Richardson, ^221111 which permits the generation 
of exact energies for systems far too large to be practicably solved by DOCI. In Fig. [^we consider a 40-site pairing 
Hamiltonian with equally-spaced levels (cp = p) at half-filling. One can see that near G « 0.2, pCCD begins to 








































FIG. 2. Differences between pCCD, pECCD, and DOCI in the dissociation of H 2 O. Left panel: Base 10 logarithms of the 
absolute value of AE (measured in Hartrees and defined in Eqn. \22\ and of 1 — 5" (defined in Eqn. |23| |. Right panel: Base 10 
logarithms of the absolute value of 1 — Sl and 1 — Sr, defined in Eqn. |27| Kinks in the pCCD results in the left panel are due 
to changes in sign. 




FIG. 3. Differences between pGGD, pEGGD, and DOGI in the dissociation of N 2 . Left panel: Base 10 logarithms of the 
absolute value of AE (measured in Hartrees and defined in Eqn. \22\ and of 1 — 5" (defined in Eqn. |23| ). Right panel: Base 10 
logarithms of the absolute value of 1 — Sl and 1 — Sr, defined in Eqn. |27| Kinks in the pGGD results in the left panel are due 
to changes in sign. 


deviate significantly from DOCI, and for G > 0.3 we can find no real solution to the pCCD equations (solutions with 
complex T amplitudes and complex correlation energies exist, but are of limited physical interest). While pECCD 
also begins to break down somewhat, it provides a much more accurate desciption of the correlations in the pairing 
Hamiltonian for large G. Note that for the half-filled forty-site Hamiltonian under consideration, a broken number 
symmetry mean field appears at Gc ~ 0.22, not coincidentally close to the value at which pCCD begins to break 
down. While pECCD is not a panacea, it provides results much superior to pCCD for the pairing Hamiltonian and 
competitive with pCCD based on the broken-symmetry mean-fielcP^ for G not too large. 

As discussed earlier, because the pECCD left-hand wave function has much better overlap with the DOCI wave 
function that does the left-hand wave function of pCCD, one would expect pECCD to yield density matrices closer to 
the DOCI density matrices than does pCCD. This is indeed the case. Figure shows fractional errors in the entries 
of the two-particle density matrix contributions {P^ Pq) for N 2 in the recoupling regime {R^-n ~ 2.2 A). Explicitly, 
we are plotting the fractional error 1 — T/Tdoci in the matrix elements of the two-particle density matrix for all 
elements Edoci larger than 10“^, since very small elements of Tdoci are unlikely to have significant contributions to 
expectation values. We see that typically, the pCCD values differ from those of DOCI by ^ 1%, while the errors in the 
pECCD density matrix elements are an order of magnitude or more smaller. Note that we have Hermitized the pCCD 
and pECCD density matrices to simplify the comparison; this has a much larger effect on pCCD than on pECCD 
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FIG. 4. Differences between pCCD, pECCD, and DOCI in the dissociation of H 2 O, pairing canonical Hartree-Fock orbitals. 
Left panel: Base 10 logarithms of the absolute value of AE (measured in Hartrees and defined in Eqn. 221 and of 1 — S' (defined 
in Eqn. 231. Right panel: Base 10 logarithms of the absolute value of 1 — S_l and 1 — Sr, defined in Eqn. 27 Kinks in the 


pCCD results in the left panel are due to changes in sign. 




FIG. 5. Energies in the 40-site pairing Hamiltonian at half-filling. Left panel: total energies. Right panel: fraction of the 
DOGI correlation energy recovered by pGCD, pEGCD, and by pCGD based on a number-broken BGS reference. Near the 
value marked Gc, number symmetry is spontaneously broken at the mean-field level, as can be seen from the BGS-based pCGD 
results. 


because the pECCD density matrix is more nearly Hermitian. For example, the Frobenius norms •\/Tr(r F^) of the 
antisymmetric parts of two-particle density matrices discussed in Fig. HI are 5.15 x IQ-^ for pCCD and 7.68 x IQ-^ 


for pFCCD. 


IV. CONCLUSIONS 

In many cases, a configuration interaction restricted to paired excitations but not restricted by excitation level 
can provide an accurate accounting for static correlation effects, particularly once the orbitals used to define the 
pairing and the reference determinant are optimized. The idea that this doubly occupied configuration interaction 
could describe many forms of static correlation is not a new one; indeed, DOCI itself was introduced over forty years 
But while pair-excited configuration interaction attracted a great deal of early interest, the method was 
basically abandoned due to its exponential scaling with system size, which renders DOCI unsuitable for practical 
calculations. 

With pair coupled cluster doubles and pair extended coupled cluster doubles, we now possess two models in the 
coupled cluster family which generally provide results almost indistinguishable from those of DOCI but with mean-field 
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loS( I Tooci') 


FIG. 6. Base 10 logarithms of the fractional differences between pCCD, pECCD, and DOCI two-particle density matrices. 
Here, dF means the difference between entries of the pCCD or pECCD density matrix and the DOCI density matrix. We have 
omitted elements of the DOCI density matrix smaller in magnitude than 10“^. 


computational cost. While the pCCD energy and right-hand wave function typically reproduce DOCI very well, the 
left-hand wave function, being a linear combination of the reference and doubly-excited determinants, is not always 
of particularly good quality, and the coincidence between pCCD and DOCI breaks down for the attractive pairing 
interaction (though note that the coefficient of the pairing-type interaction in the physical electronic Hamiltonian is 
positive, and can only become negative upon renormalization with some Hamiltonian transformation). Both these 
difficulties can be ameliorated by using pair extended coupled cluster doubles, which seems to offer much more accurate 
results for attractive interactions and a much superior left-hand wave function (and therefore superior density matrices 
and expectation values). The computational scaling of pECCD is the same as pCCD, namely 0{N^), though the 
pECCD energy and amplitude equations are significantly more complicated and the method is accordingly rather more 
expensive than is pCCD itself. Nonetheless, the mean-field computational scaling permits routine pECCD calculations 
on systems with hundreds of basis functions, for which we can reliably anticipate getting results of essentially DOCI 
quality for both the energy and for other observables. We hope, then, that pECCD will be a valuable tool for the 
description of strongly correlated systems, particularly when pairing interactions become attractive and G is not too 
large, and that it will form an excellent starting point from which to add correlations which break pairs. 

There are, of course, important drawbacks of pECCD as well. The most significant is that pECCD does not 
include dynamic correlation. One can attempt to fix this with the addition of a simpl e de nsity functional correlation 
energy,!^ or by relaxing the restriction to seniority zero, whether by freezing amplituded^ or, potentially, by including 
a perturbative account for higher seniority sectors. Moreover, the restriction to seniority zero is sometimes too severe; 
some systems simply require higher seniority sectors for the description of strong correlation effects, as is seen, for 
example, in the dissociation of N 2 .^ We may hope to treat these problems by lifting the restriction that the pairing 
scheme pairs the two spinorbitals corresponding to the same spatial orbital. This might even open the door to the 
use of broken spin symmetry in seniority zero methods for strongly correlated systems. 
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Appendix A: Useful Derivations 


Here, we sketch a few derivations which the reader may find helpful but which are not necessary for understanding 
the thrust of the manuscript. 
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1. Derivation of ° 


Recall that is the portion of the Hamiltonian which preserves seniorities of individual spatial orbitals. Here, 

we wish to derive the expression quoted in Eqn. |12[ 

Let us start, then, with the one-electron part of the Hamiltonian. In order to preserve orbital seniority, when we 
remove an electron from orbital p we must return it to orbital p. Thus, we have 



^'^{P\h\p)'^al^ Qp^ 

(Ala) 

pq a 

p <7 



= ^{p\h\p) X’^P- 

p (7 

(Alb) 

where the spinorbital number operators are 

= ^^\h\p)Np 

p 

(Ale) 

riro = 
per 


(A2) 


and their sum is the spatial orbital number operator Np given in Eqn. 14a 

The two-electron part of the Hamiltonian is slightly more complicated. We could remove two electrons from orbital 
p and place them in orbital q, or we could remove one electron from orbital p and another from orbital q ^ p, in 
which case we must place electrons back in orbitals p and q; in the latter case, we must include the possibility of an 
exchange where the first electron is removed from p but placed in q. All told, we have 


\ XI al^ al^ X 4. 4, 

pqrs (TT) pq err} 

X '^iP<i\v\pq) “p- 

p^q at] 

“Pr, 

p^q aq 

=\ ^{ppHqq) X“L “L “-Jij “9- 

pq cr 

+\ '^{pq\'v\pq) X"p- 

p^q a-q 

-\Y.^P<iM^p) X“L “p. 

p^q (jq 

= ^{pp\v\qq) Pi Pq 

pq 

+l '^{pq\v\pq) NpNq 

Py^q 

[X ^P- ^9- +Sp Sq + Sp S+ 
p^q \ a 


(A3a) 


(A3b) 


(A3c) 


where the spin index tr is the opposite of the index a and where we have made use of the pair creation operator Pi 
given in Eqn. [^and the adjoint operator Pq and have introduced the spin raising and lowering operators 


5+ = al^ ap^ 

(A4a) 

Q — t 

^p ~ ®P4. ®Pt ■ 

(A4b) 

^pt “ 2 “^P ’ 

(A5a) 

rip^ = 2 P^p ~ ^p 

(A5b) 


To simplify further we note that 



where 


— 2 (®Pt ®P4. “ 2 

^1/ 

and <T^ is a Pauli matrix. This means that 

(7 

whence 

^ E “L 4r, ^ H ^9 + ^ E (P'i'kbg) - (wkkp)) Np Nq 

pqrs ( 77 ] pq p^q 

- I (2 ^P ^P + “^P '^ 9 ” + 5'p -5'+) 

P7^9 

“ E/ ^P9 -fp ^9 + ^ E/ ^P9 -^P -^9 + E/ ^^^9 Sp ' Sq, 

pq p^q p^q 

in terms of the spin vector operators defined in Eqn. |14b| so that overall, we may make the replacement 

H ^ E hp Np + ^ ^ ^ Wpq Np Nq + 'y ^ Vpq Pp Pq + E Kpq Sp ■ Sq, 
p p¥=q PI p¥^q 

as desired. 


(A6) 


(A7) 


(A8a) 


(A8b) 


(A9) 


2. Density Matrices 

In order to compute the pECCD energy, we need to evaluate density matrix elements. In other words, we need to 
take the expectation value 


(77o'''=°)pECCd = (0|e^e-^i7o^^=°e^|0), (AIO) 

where we recall that |0) is the single-determinant reference. As noted in Ref. [HI the density matrices of seniority zero 
methods are sparse, and in fact the only non-zero entries of the full one- and two-particle density matrices can be 
constructed from (Ap)pECCD, {NpNq)p^cCD, and P 5 )pECCD- 

Computing the density matrix elements is simplified by constructing similarity transformations of the number 
operators Np and the pair creation and annihilation operators P^ and Pq. We find 


e-^ A^a = A^a + 2 E p1 

(Alla) 

e-^ N,e^ = N,-2 Y,ka Pi 

(Allb) 

e-^Ple^ = Pl, 

(Allc) 

e-^ Ple^ = Pl+J2 Ub Pi [N, - 1) - E Pi Pi 

(Alld) 

b a^b 


e“^ e^ = Pa + E “ E ^9“ pI P^ Pp 

(Alle) 

i i=jtj 


e-^P,e^ = P„ 

(Allf) 


where we have taken advantage of the commutator expansion 
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and have used the SU(2) commutation relationships given in Eqn. 15 and transcribed here for convenience: 


[PpiPl] — ^pq (1 ^p) ) 

[Np^Pq] = -‘^SpqPq. 


(Al3a) 

(A13b) 


We have added a few summation index restrictions which are unnecessary (but permitted) in view of the nilpotency 
of the pair creation and annihilation operators, simply because they may help cla rify the origins of the summation 
restrictions and factors such as Sab and Sij in the hnal density matrices of Eqn. 18 
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